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1 Introduction 



There has been recent interest, particularly in the systems biology literature, in the study of sym- 
metry invariances of responses of dynamical systems. The paper [1] obtained sufficient characteriza- 
tions of symmetry invariance using a notion of equivariance, and this characterization was shown to 
be necessary as well as sufficient in [2] . Both [I] and [2 J sketched how to extend the results to motile 
systems that explore space, so long as the "motor dynamics" depends only on an invariant response. 
Specifically, these results predicted that E. coli bacteria would produce scale-invariant searches, 
meaning that distributions of bacteria, even under non-uniform and time-varying chemoeffector 
fields, should be invariant under any rescaling of the input field. This prediction was subsequently 
experimentally verified in [3j. In this note, we remark that, for a velocity-jump Markov model, the 
PDE for the evolution of densities (or normalized concentrations) in time inherits the symmetry- 
invariance property from individual behaviors. Although not at all surprising, this provides further 
theoretical justification for passing from individual-based models to population predictions. 



2 Symmetries and equivariances 

We review the general setup in [TJ [2] . Consider dynamical systems with inputs and outputs [I] , 

x = f(x,u), y = h(x,u). (1) 

The functions /, h describe respectively the dynamics and the read-out mapQ Equation is 
shorthand for 

dx 

— (t) = f(x(t), u(t)) , y(t) = h(x(t),u(t)) ■ 

Here, u = u(t) is a generally time-dependent input (stimulus, excitation) function, x(t) is an n- 
dimensional vector of state variables, and y{t) is the output (response, reporter) variable. States, 
inputs, and outputs are constrained to lie in particular subsets X, U, and Y respectively, of Euclidean 
spaces R n ,K m ,K 9 . 

We assume that for each piecewise-continuous input u : [0, oo) — > U, and each initial state £ G X, 
there is a unique solution x : [0,oo) — > X of ([I]) with initial condition x(0) = £, which we write 

"The results in [5] were stated for h not directly dependent on u, but the theory is the same in the more general 
case of u-dependence, as was also remarked there. 



as ip(t, £, u), and we denote the corresponding output y : [0, oo) — > Y, given by h(ip(t,£,u),u(t)), 
as ip(t,£,u). We also assume that for each constant input u(t) = u, there is a unique solution 
x = o~(u) of the algebraic equation f(x,u) = 0. Often one also assumes that this steady state is 
globally asymptotically stable (GAS): it is Lyapunov stable and globally attracting for the system 
when the input is u(t) = u: limt_»oo ip(t, £, u) = o~(u) for every initial condition ( £ X. The GAS 
property is not required for the results to follow, however. 

If X is an open set, or the closure of an open set, in IR n , the system is said to be analytic if / 
and h are real-analytic (can be expanded into locally convergent power series around each point) 
with respect to x, and irreducible if it is accessible and observable. 

An accessible system is one for which the accessibility rank condition holds: J-la(xo) = K n for 
every xq G X, where J-la is the accessibility Lie algebra of the system. Intuitively, this means that 
no conservation laws restrict motions to proper submanifolds. For analytic systems, accessibility is 
equivalent to the property that the set of points reachable from any given state x has a nonempty 
interior; see a proof and more details in the textbook [4]. An observable system is one for which 
ip(t, xq, u) = ip(t, xq, u) for all u, t implies xq = x~q. Intuitively, observability means that no pairs of 
distinct states can give rise to an identical temporal response to all possible inputs. For analytic 
input-affine systems, observability is equivalent to the property that any distinct two states can 
be separated by the observation space; see [I], Remark 6.4.2 for a proof and discussion. In the 
context of applications to biomolecular systems, analyticity and irreducibility are weak techinal 
assumptions, often satisfied. 

Adaptation, invariance, and equivariances 

Definition 2.1 The system ([!]) perfectly adapts to constant inputs provided that the steady-state 
output h(a(u),u) equals some fixed yo G Y, independently of the particular input value u G U. □ 

That is, the steady-state output value is independent of the actual value of the input, provided 
that the input is a constant (a step function). 

Invariance will be defined relative to a set V of continuous and onto input transformations tt : U — > 
U. For each input u(t) and tt £ V, we abuse notation and denote by "iru" (even if tt is nonlinear) 
the function of time that equals ir(u(t)) at time t. (The continuity assumption is only made in 
order to ensure that ttu is a piecewise continuous function of time if u is. The ontoness assumption, 
that is, 7rU = U, and can be weakened considerably: it is only used in in the main theorem in order 
to prove that a system x = f(x, ttu), y = h(x, ttu) is irreducible if the original system is irreducible, 
but far less than ontoness is usually required for that.) 

An example is scale invariance, in which U = M>o and V = {u i— > pu,p > 0}. (Scale invariance is 
sometimes called "fold-change detection" (FCD), since the only changes that can be detected in a 
response are those due to different fold-changes in inputs.) 

Definition 2.2 The system has response invariance to symmetries in V or, for short, is V- 
invariant if 



Under the assumption that the action of V is transitive, i.e., for any two u, v G U, there is some it 
such that v = ttu, "P-invariance implies perfect adaptation, because the outputs in ^ must coincide 
at time zero, and any two inputs can be mapped to each other. 





holds for all t > 0, all inputs u = u(t), all constants u, and all transformations ir £ V. 



□ 
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Definition 2.3 Given a system ([!]) and a set of input transformations V, a parametrized set of 
differentiable mappings {p w : X — > Xj^-p is a V -equivariance family provided that, for each n: 

f(p 7r (x),iru) = (p 7r )*(x)f(x,u) and h(p n (x), nu) = h(x, u) (3) 

for all x € X and u G U, where (p w )* denotes the Jacobian matrix of p n . If ^ holds, the system is 
said to be pjr-equivariant under the input transformation ir. □ 

The first part of Equation ^ is a first order quasilinear partial differential equation on the n 
components of the vector function p n , for each u G U, and one may solve such equations, in 
principle, using the method of characteristics. The second part of Equation ^ is an additional 
algebraic constraint on these components. Observe that the verification of equivariance does not 
require the computation of solutions ip(t, a(7ru),7Tu). We omit the subscript tt when clear from the 
context. 

The main result in [2] is as follows. 

Theorem 1 An analytic and irreducible system is V -invariant if and only if there exists a V- 
equivariance family. I 

Remark 2.4 An interesting consequence of this theorem is that, if "P-invariance holds, then a 
stronger property holds as well, namely that 

Tp(t,X,u) = 1p(t, p(x), TTU) 

is valid for all t > 0, all inputs u, all transformations tt £ V, and every initial state x (not necessarily 
a steady state). Another interesting fact, which follows from the proof of the theorem, is as follows. 
Suppose that we define a "weakly invariant" system as one for which there exists some constant u 
such that ([2]) holds: ip(t,a(u),u) = ip(t,a(iru),Tru) for all inputs u and all i > (instead of asking 
that this holds for every u). Then, "weak invariance" implies the existence of an equivariance, and 
hence also invariance. The irreducibility property plays a subtle role in these facts. □ 



3 Symmetry-invariant steering 

We consider next a motile vehicle or organism which explores a space while measuring the "inten- 
sity" of an input cue (such as a chemoeffector or light). The sensed input at time t and position r 
is U(t, r), where r = r(t) is the current position of the vehicle. The current position r{t) is derived 
from the output y(t) of a system ([!]), through a computation that takes into account the dynamics 
of the motor and steering mechanisms. 

Deterministic models for such mechanisms are sometimes appropriate, and one was described in [TJ 
12]. An easy argument for that deterministic model shows that, if y is invariant under symmetries 
in inputs, then positions r(t) will be invariant under symmetry transformations on U. 

It is often the case that a more accurate description is one in which the output y(t) drives a 
stochastic, not a deterministic, steering mechanism: the subsystem producing the location r(t) is 
subject to randomness. 

An important instance of this is bacterial E. coli chemotaxis, where y(t) represents a signal, the 
level of phosphorylated protein CheY, which serves to bias the random switches between tumbling 
and swimming ("run") modes. Specifically, let us consider the Tu-Shimizu-Berg E. coli chemotaxis 
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model [5j, which may be formulated, for realistic parameters and input levels, as follows: m = Fo(y), 
y = h(m, u) = G(u/e am ), where Fq is a decreasing function which crosses zero at some value y = yo 
(and G is a suitable function whose precise form is immaterial for establishing symmetry). Letting 
x ■= e am and F = aF, we may transform this system into a "nonlinear integral feedback" form, 

x = xF(h(x,u)) 
h(x, u) = G (u/x) . 

For this system, homogeneity of f(x,u) = xF(h(x,u)) implies scale invariance, since the unique 
solution of the equivariance PDE is p(x) = px, for the scaling symmetry u h-> pu. Based on this veri- 
fication of scale-invariant behavior, pQ predicted the invariance of distributions of bacteria locations 
under scalings of chemoattractant fields. This prediction was subsequently verified experimentally 
in [3] by means of molecular level analysis of intracellular signaling (FRET experiments) as well as 
measurements of swimming behavior at the level of individual cells and populations (in microfluidic 
environments) . 

A simple numerical simulation serves to illustrate the point. This simulation uses (with no change 
in parameters) the SPECS agent-based model for E. coli chemotaxis that was developed in [6jj^] In 
this simulation, cells are allowed to swim in a rectangular channel that is 2000^m long and 400/im 
wide, and data is collected in bins of size 20 (so, there are 100 bins along the long axis). The ligand 
gradient is stationary and linear (see below for boundary values) along the length and constant along 
the width. We simulated 1000 cells, all initially placed at the middle (at length 1000, i.e. bin 50), 
and plotted the marginal distributions (along the long axis on which the chemoattractant varies). 
Since there behavior is random, the averages of several (five) trials under each of the conditions are 
shown. These average histograms are plotted for the cell distribution at time t = 500. The blue and 
green histograms in Fig. [T] represent, respectively, results for cells pre-adapted to a concentration 
250 (units are /jlM), and linear gradient 200 . . . 300, and cells pre-adapted to a concentration 375, 
and linear gradient 300 . . . 450 (a scale change by p = 1.5). As expected, the distributions are very 
similar. As a control, we also plotted the results of using, once again, a linear gradient 300 . . . 450, 
but now pre-adapting cells to a concentration of 250. Since the initial state is not matched, there 
is no reason for invariance. Indeed, the resulting red histogram is very different from the previous 
ones. 
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Figure 1: Simulations using SPECS code 
We thank Y. Tu for making this code available. 
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One may mathematically formalize probabilistic behavior, and show symmetry-invariance of search 
under randomness, in several possible ways. For instance, in [2] a simple result was presented on 
symmetry-invariance search based on pathwise equality of stochastic processes. We describe next a 
different approach, that employs the formalism of velocity-jump processes [7] with added internal 
dynamics [8]. 

We wish to model motions in a space M. N (typically N = 1,2,3; and we assume for simplicity 
that motion can occur on the entire space) of individuals (bacteria, vehicles, etc) whose internal 
dynamics are described by the states x in and which change velocities as a function of the 
output y. To avoid confusion with the variable x used for the internal state, we use the letter "s" 
to denote points in the space ~R N in which movement occurs. The input u = u(t, s) represents an 
external signal present at time t in location s. The subset V C M. N denotes the space of possible 
velocities. 

We assume that the system can instantaneously change orientations. (For E. coli bacteria this 
would mean that we are ignoring tumble durations.) 

The concentration at time t of individuals present at time t in location s and having internal 
state x and velocity v is denoted by c(t, s, v, x). We interpret c(t, s, v, x) ds dv dx as the number of 
individuals located between s and s + ds, having velocity between v and v + dv, and whose internal 
state is between x and x + dx. Normalized by the total number of individuals, one may also think 
of probability density, at each time t. 

We assume that velocities change at random. The times at which velocities jump are controlled by 
a Poisson process with intensity X(y). Given that a jump in velocity occurs, which particular new 
velocity is picked is itself the result of a random choice; the kernel T y (v, v' , y) gives the probability 
of a change in velocity from v' to v. Since T is a probability density, j v T y (v,v') dv = 1 for every 
y. Notice that, just as with the jump instants, the kernel also depends on the state only through 
the output y. 

Then the evolution (transport, Fokker- Planck, or forward Kolmogorov) equation for c = c(t, s, v, x) 
is: 

dc f 

— + V s -cv + V x -cf = -\{y)c + J^X(y)T y {v,v')c(t,s,v',x)dv' (4) 

where most arguments have been omitted for simplicity, but understood as holding for all (t, s, v, x). 
(More generally, the right-hand side could be replaced by a more complicated discrete rate of change, 
if jumps are governed by a non-Poisson process.) The input at location s and time t is U(t, s), and 
it appears in these equations through the vector field / in Q. 

Sometimes it is useful to view Q as a set of partial differential equations indexed and coupled by 
the velocities v. For example, when N = 1 and there is a constant speed vq > 0, V = {— vq,vo} 
is a two-element set which provides the orientation of movement, T y (v,v') = 1 (there is only one 
possible jump, namely a reversal of direction), and Q describes a telegraph-type process: denoting 
c + (t, s, x) = c(t, s, vq, x) and c~(t, s, x) = c(t, s, —vo, x), Q can then be thought of as set of coupled 
partial differential equations, one for dc + /dt and one for dc~ /dt: 

~dT + V0 ^x~ + Vx ' fc+ = A (y)[- c+ + c_ ] 
~m~ v °~dx~ + Vx ' fc ~ = A (y)[c + -c - ]. 

The reference [8] discusses mathematical aspects of the PDE Q, which will not be discussed here. 
We focus, purely formally, on symmetry invariance. 
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Let us assume given tt and an associated equivariance p = p n , so that ((3j) holds: 

f(p(x),iru) = p*(x)f(x,u) and h(p(x), iru) = h(x, u) 

for all x E X and u E U, where p* denotes the Jacobian matrix of p. We will also make the following 
assumption on the divergence of /: 

(V x - f)(p(x),7T U ) = (V x -f)(x,u) (5) 

for all x E X and m£D. This property is automatically satisfied for most of the examples treated 
in [2], since in these examples, which are for scale invariance ttu = pu, p is a linear mapping. In 
general, if p(x) = Rx for a matrix R, then the equivariance condition f(Rx, ttu) = Rf(x, u) implies, 
taking Jacobians, that f*(Rx, ttu) = Rf*(x, u)R~ l . Since two similar matrices have the same trace, 
and V x • / is the trace of the Jacobian of /, it follows that ^ is valid. 

Our main observation is that the same distribution of individuals will result if the input field U is 
replaced by nil, provided that the internal states are transformed by p. A precise statement is as 
follows. 

Theorem 2 Suppose that c satisfies with respect to an input field U. Define 

c(t, s, v, x) = c(t, s, v, p~ 1 (x)) . 
Then c satisfies with respect to the input field irU . 

Proof. We start by writing all the arguments in Q explicitly: 

Oc 

— (t, S, V, x) + (V a • Ti)(t, S, V, x) + (V x • r 2 )(t, S, V, x) 

= —X(h(x,U(t,s)))c(t,s,v,x) 

+ I X(h(x, U (t, s)))T h ( x>u{ttS)) (v, v')c(t, s, v', x) dv' 

where 



v 



r%(t,s,v,x) = c(t,s,v,x)v 
T 2 (t,s,v,x) = c(t,s,v,x)f(x,U(t,s)) . 

Since this equation must hold for all x, it holds also when p _1 (x) is replaced for x, in other words 
it is also true that 

— (t, s, v, p-^x)) + (V s • ri)(t, s, v, p-^x)) + (V x ■ r 2 )(t, s, v, p-^x)) 
= -\{h{p-\x),U{t, s)))c(t, s, v, p-\x)) 

+ / KKp~ l { x ),U{t,s)))T h{p -i {x))U{t)S)) (v,v')c{t,s,v\ p~ l {x)) dv' 

J \f 

for all t, s, v, x. From the definition of c and the property h(p(x),iru) = h(x, u), which implies that 
h(x,iru) = h(p~ 1 (x),u) for all u, we conclude that: 

— {t, s, v, x) + (V s ■ ri)(t, s, v, p-^x)) + {V x ■ r 2 )(t, s, v, p-^x)) 
= —X(h(x,irU(t,s)))c(t,s,v,x) 

+ / \(h(x,irU(t,s)))T h{x!7TU{ttS)) (v,v')c(t,s,v',x) dv' . 
Jv 
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It will follow that c is a solution of Q with respect to the input field ttII provided that we show: 

(V x -r 2 )(t,s,v, /0 - 1 (x)) = (V x -f 2 ){t,s,v,x), 

where 

Ti(t, s,v,x) = c(t,s,v,x)v 
F2(t,s,v,x) = c(t,s,v,x)f(x,irU(t,s)). 

Since c(t, s,v,x) = c(t,s,v, p (x)), the equality for V s is clear. We are left to show the equal- 
ity for V x . We have, fixing t, s, v and letting F(x) = f(x,U(t,s)), C(x) = c(t, s,v,x), F(x) = 
f(x, irU(t, s)), and C(x) = c(t, s, v, x) = C(p~ 1 (x)): 

V x ■ [CF](x) = (dC/dx)(x)F(x) + C(x)(y x ■ F){x) 

= (dC/dx)(p- l {x)){p- l )*(x)F(x) + C(p^(x))(V x ■ F)(x) 
= (dC/dxXp-Hx^p-^x))]- 1 ^) + C(p-\x))(V x ■ F)(x) 
= (dC/dx)(p~\x))F{p'\x)) + C{p-\x))(V x ■ F)(p-\x)) 
= V x -[CF](p- 1 (x)), 

where we have used that f(p(x),Tru) = p*(x)f(x,u), and thus also F(x) = p*(p~ 1 (x))F(p~ 1 (x)), 
as well as the divergence property ([5]). I 

In applications, one is often interested in the distribution of positions irrespective of internal states 
x and velocities v. 

Q(t,s)= / / c(t, s,v ,x) dpx(x) dpv(v) 
Jx Jv 

where px an d pv denote appropriate measures on X and V (and we assume that c is integrable). 
Take the density corresponding to nil, c(t, s, v, x) = c(t, s, v, / o~ 1 (x)), and its marginal 

Q(t,s)= / / c(t,s,v,x)dpx(x)dp v {v) . 
Jx Jv 

This is the same as f x f v c(t, s,v,x)r(x) dpx(x) dpv{v), where r(x) = 1/ detp*(x). In the special 
(but usual in examples) case that p is linear, r is a constant, so Q(t, s) = rQ(t, s). It follows that 
the normalized densities are equal: 

Q(t,s) = Q{t,s) 
fQ(t,(r)da ' ' fQ(t,a)da ' 

Alternatively, one could introduce a new measure dpx(x) = r(x)px, and define Q using this new 
measure, for all times t and space positions s, so that Q(t, s) = Q(t, s). 
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